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A formalism is introduced for the non-perturbative, purely numerical, solution of the reduced 
Rayleigh equation for the scattering of light from two-dimensional penetrable rough surfaces. As 
an example, we apply this formalism to study the scattering of p- or s-polarized light from two- 
dimensional dielectric or metallic randomly rough surfaces by calculating the full angular distribution 
of the CO- and cross-polarized intensity of the scattered light. In particular, we present calculations 
of the mean differential reflection coefficient for glass and silver surfaces characterized by (isotropic 
or anisotropic) Gaussian and cylindrical power spectra. The proposed method is found, within the 
validity of the Rayleigh hypothesis, to give reliable results. For a non-absorbing metal surface the 
conservation of energy was explicitly checked, and found to be satisfied to within 0.03%, or better, 
for the parameters assumed. This testifies to the accuracy of the approach and a satisfactory 
discretization. 

PACS numbers: 42.25.-p, 41.20.-q, 78.20.-e, 78.20.Bh 



I. INTRODUCTION 

Wave scattering from rough surfaces is an old disci- 
pline which keeps attracting a great deal of attention 
from the scientific and technological community. Several 
important technologies in our society rely on such knowl- 
edge, with radar being a prime example. In the past, the 
interaction of light with rough surfaces was often con- 
sidered an extra complication that had to be taken into 
account in order to properly interpret or invert scattering 
data. However, with the advent of nanotechnology, rough 
structures can be used to design novel materials with tai- 
lored optical properties. Examples include: metamateri- 
als [HE], photonic crystals [3], spoof plasmons optical 
cloaking [S}^ , and designer surfaces |8i i9j . These devel- 
opments have made it even more important to have avail- 
able efHcient and accurate simulation tools to calculate 
both the far- and near-field behavior of the scattered and 
transmitted fields for any frequency of the incident radi- 
ation, including potential resonance frequencies of the 
structure. 

Lord Rayleigh was the first to perform systematic stud- 
ies of wave scattering from rough surfaces when, in the 
late 1800s, he studied the intensity distribution of a wave 
scattered from a sinusoidal surface [TOl [11] . More than 
three decades later, Mandel'shtam studied light scatter- 
ing from randomly rough surfaces [TJ thereby initiat- 
ing the field of wave scattering from surface disordered 
systems. Since the initial publication of these seminal 
works, numerous studies on wave scattering from ran- 
domly rough surfaces have appeared in the literature [T51 - 
119] , and several new multiple scattering phenomena have 
been predicted and confirmed experimentally. These 
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phenomena include the enhanced backscattering and en- 
hanced transmission phenomena, the satellite peak phe- 
nomenon, and coherent effects in the intensity-intensity 
correlation functions [1^24) . 

These studies, and the methods they use, can be cat- 
egorized as either perturbative or purely numerical (and 
non-pertubative). While the former group of methods 
is mainly limited to weakly rough surfaces, and there- 
fore have limited applicability, the latter group of meth- 
ods can be applied to a wider class of surface rough- 
nesses. Rigorous numerical methods can in principle be 
used to study the wave scattering from surfaces of any 
degree of surface roughness. Such simulations are rou- 
tinely performed for systems where the interface has a 
one-dimensional roughness, i.e., where the surface struc- 
ture is constant along one of the two directions of the 
mean plane |19l I25j . However, for the practically more 
relevant situation of a two-dimensional rough surface, the 
purely numerical and rigorous methods are presently less 
used due to their computationally intensive nature. The 
reason for this complexity is the fact that for a randomly 
rough surface there is no symmetry or periodicity in the 
surface structure that can be used to effectively reduce 
the simulation domain. For a periodic surface, it is suf- 
ficient to simulate a single unit cell, while for a random 
surface the unit cell is in principle infinite. 

A wide range of simulation methods are currently avail- 
able for simulating the interaction of light with mat- 
ter, including the finite-difference time-domain (FDTD) 
method the finite-element method (FEM) [?fl[^ . 
the related surface integral equation techniques also 
known as the boundary element method (BEM) or 
the method of moments (MoM) |29ll33) . the reduced 
Rayleigh equation (RRE) technique [HI [MSO], and 
spectral methods [5^ . 

The FDTD and FEM methods discretize the whole 
volume of the simulation domain. Due to the complex 
and irregular shape of a (randomly) rough surface, it is 
often more convenient, and may give more accurate re- 
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suits (for the same level of numerical complexity) [H] , to 
base numerical simulations on methods where only the 
surface itself needs to be discretized. This is the case, 
for example, for the surface integral technique and the 
reduced Rayleigh equation methods. 

The reduced Rayleigh equation is an integral equation 
where the unknown is either the scattering amplitude 
or the transmission amplitude. In the former (latter) 
case, one talks about the reduced Rayleigh equation for 
reflection (transmission) . For reflection this equation was 
originally derived by Brown et al. |34j . and subsequently 
by Soubret et al. [39l HQ] . Later it has also been derived 
for transmission [32] and film geometries [351 1131 SI] ■ 

In the past, the surface integral technique has been 
used to study light scattering from two-dimensional ran- 
domly rough, perfectly conducting or penetrable sur- 
faces [33 [33 US]- However, to date, a direct numeri- 
cal and non-perturbative solution of the two-dimensional 
reduced Rayleigh equation has not appeared in the lit- 
erature, even if its one-dimensional analog has been 
solved numerically and has been used to study the scat- 
tering from, and transmission through, one-dimensional 
rough surfaces [35H37] . The lesson learned from the one- 
dimensional scattering studies reported in Refs. [35H57] 
is that simulations based on a direct numerical solution 
of the reduced Rayleigh equation may give accurate non- 
perturbative results for systems where alternative meth- 
ods struggle to give the same level of accuracy. Moreover, 
the reduced Rayleigh method also requires less memory 
for the same surface dimensions when compared to, e.g., 
the rigorous surface integral technique. 

The main aim of this paper is to present a numeri- 
cal method and formalism for the solution of the two- 
dimensional reduced Rayleigh equation for reflection. 
While we exclusively consider reflection, the formalism 
for transmission will be almost identical, and the result- 
ing equation will have a similar form as for reflection. 
Additionally, the equation for transmission or reflection 
for a film geometry, i.e., for a film of finite thickness on 
top of a substrate, where only one interface is rough, will 
also have a similar form. The method presented will be 
illustrated by applying it to the study of the scattering of 
p- or s-polarized light from two-dimensional metallic or 
dielectric media separated from vacuum by an isotropic 
or anisotropic randomly rough surface. 

This paper is organized as follows: First, in Sec. [Il] we 
present the scattering geometry to be considered. We will 
then present some relevant scattering theory, including 
the reduced Rayleigh equation for the geometry under 
study (Sec. Ill), followed by a detailed description of how 
the equation can be solved numerically (Sec. IV). Next, 



we will present some simulation results obtained by the 
introduced method (Sec.|v]). We then discuss some of the 
computational challenges of this method (Sec. VI), and. 




FIG. 1. (Color online) A sketch of the scattering geometry as- 
sumed in this work. The figure also shows the coordinate sys- 
tem used, angles of incidence (6o,(j)o) and scattering {9s,(j)s), 
and the corresponding lateral wavevectors ky and qy , respec- 
tively. 



II. SCATTERING GEOMETRY 

We consider a system where a rough surface sepa- 
rates two regions. Region 1 is assumed to be vacuum 
(ei = 1), and region 2 is filled with a metal or dielec- 
tric characterized by a complex dielectric function S2{ttj), 
where the angular frequency is a; = 27rc/A, with A being 
the wavelength of the incident light in vacuum and c the 
speed of light in vacuum. The height of the surface mea- 
sured in the positive direction from the xia::2-plane is 
given by the single- valued function x^ = C(x||), where 
X|[ = (a;i,a;2,0), which is assumed to be at least once 
differentiable with respect to Xi and X2- Angles of inci- 
dence (6*0, (po) and scattering {9s, 4>s) are defined positive 
according to the convention given in Fig. [ij 

In principle, the theory to be presented in Sec. |III| can 
be used to calculate the scattering of light from any sur- 
face, provided it is not too rough. However, in this paper, 
we will consider randomly rough surfaces where C(x||) 
constitutes a stationary random process defined by 



(C(x||))=0, 
(C(x||)C(x||'))='52l^(x|| 



(1) 



where the angle brackets denote an average over an en- 
samble of surface realizations. In writing Eqs. ([!]) we 
have defined the root-mean-square height of the surface, 

1/2 

5 = (C^(x||)) , and W{yi\\ — X||') denotes the height- 
height auto-correlation function of the surface, normal- 
ized so that W^(0) = 1 [19 . According to the Wiener- 
Khinchin theorem [47j , the power spectrum of the surface 
profile function is given by 



finally, in Sec. VII we draw some conclusions. 



5(ki| 



d^a 



Vl^(x|[ ) exp (— ikii • xii 



(2) 
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The power spectra that will be considered in this work 
are of either the Gaussian form 1331 



.9(k|| 



'7raia2 exp 



kiaf 



(3) 



where {i — 1,2) denotes the lateral correlation length 
for direction i, or the cylindrical form [35] 



5(^11 ) = 



47r 



fc2 



(4) 



where fc|| = |k|||, 9 denotes the Heaviside unit step 
function, and k± are wavenumber cutoff parameters, 
with fc„ < fc+. The cylindrical form in Eq. Q is a 
two-dimensional generalization of the power spectrum 
used in the experiments where West and O'Donnell con- 
firmed the existence of the enhanced backscattering phe- 
nomenon for weakly rough surfaces |21j . 



III. SCATTERING THEORY 

We consider a linearly p- or s-polarized plane wave 
which is incident on the surface from region 1, with the 
electric field given by E'"'=(x;t) = E'"=(x|w) exp(-iwt) 
where 



E"-(x|u;) = 5'"^(k||)exp [ik,i • xy - iai(fc,i)x3] 
with 



5'""(k||) 



k||ai(fc||) + xafcii 
X3xk||)£r(k||), 



(5a) 



(5b) 



and 



ai(fc||) = 



-fcf?, Reai > 0, Imai > 0. (5c) 



Here, and in the rest of the paper, a caret over a vector 
indicates a unit vector. The expression s in front of the 
amplitudes f^"'^(kj|) (a = p, s) in Eq. (5b) correspond 
to unit polarization vectors for incident light of linear 
polarization a. Moreover, ky = (fci,fc2,0) denotes the 
lateral component of the wave vector k = k|| — a(fc||)x3. 
When the lateral wavenumber satisfies fcy < w/c, as will 
be assumed here, ky is related to the angles of incidence 
according to 



kii ~ — sin ^0 (cos ^q, sin 0) 
c 



where c denotes the speed of light in vacuum and 
and 00 9'i'6 the polar and azimuthal angles of incidence, 
respectively (Fig. [T]). When writing the field of inci- 
dence, E'"'^(x; t), a time harmonic dependence of the form 
exp(— iwt) was assumed. A similar time dependence will 
be assumed for all field expressions, but not indicted ex- 
plicitly. 



Above the surface roughness region, i.e., for ^3 > 
maxC(xy), the scattered field can be written as a super- 
position of upwards propagating reflected plane waves: 



E"=(x|a;) 



^^^(qii) 



(2^)2 (7a) 
X exp [iqii • xn + \ai [qu )x^] , 



where 



r^(q||) = -[q||ai(9||)-X3q||]5;^(q||) 



(x3Xq||)fr(q||). 



(7b) 



The integration in Eq. ( |7a[ ) is over the entire plane, 
including the evanescent region gy > ui/c. Therefore, 
both propagating and evanescent modes are included in 
E"'=(x|a;). 

We will assume that a linear relationship exists be- 
tween the amplitudes of the incident and the scattered 
fields, and we write (for a = p, s) 

f-(qy)= Rc.M\\H)^'riH)- (8) 

/3=p,s 

Here we have introduced the so-called scattering am- 
plitude i?Q^(qy|ky), which describes how incident /3- 
polarized light characterized by a lateral wave vector 
ky is converted by the surface roughness into scattered 
light of polarization a and lateral wave vector qy . When 
q\\ < uj/c, the wave vector qy is related to the angles of 
scattering [Os,4>s) by 



sin 9s (cos 4>s , sin 0s , 0) . 



(9) 



Below the surface region, i.e., for 2:3 < min C(x||), the 
transmitted electric field can be written as 



E"'(x|||t^) 



d2p| 



(27r) 
X exp fip 



2^"(P||) 



(10a) 



ia2(P||)a;3] 



with 



- [P||a2(P||) +X3Py] fp'(Pll) 



1 c 

v/£2(w)^ 

(x3Xpy)5r(pn) 



(10b) 



(6) In writing Eqs. ( 10 1 we have introduced wave vectors of 



the transmitted field p = Py — a2(?'||)x3, where 
a2{p\\ 



-pI 



(11) 



Rea2 > 0, ImQ;2 > 0. 



In complete analogy to what was done for reflection, a 
transmission amplitude TQ^(py|ky) may be defined via 
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the following linear relation between the amplitudes of 
the incident and transmitted fields (a = p, s) 



/3=P,s 



|k|t)4"'(k||). 



(12) 



Since the form of the electric fields given by Eqs. ([5]), 
([t]), and ( 10 1 apply far away from the surface region, they 



are referred to as the asymptotic forms of the electric 
field. These equations automatically satisfy the bound- 
ary conditions at infinity. 

In passing we note that once the incident field has been 
specified, the scattered and transmitted fields are fully 
specified outside the surface roughness region if the re- 
flection [i?Q^(q|| |k||)] and transmission [Tq^(p|| jky)] am- 
plitudes are known. We will now address how the reflec- 
tion amplitude can be calculated. 



A. The Rayleigh Hypothesis 

Above the surface, i.e., in the region 2:3 > maxC(x|[), 
the total electric field is equal to the sum of the inci- 
dent and the scattered field, E'"'^(x|a;) -I- E"'=(x|w). Be- 
low the surface, in the region 2:3 < minC(x|j), it equals 
the transmitted field, E*''(x|a;). In the surface rough- 
ness region, minC(x||) < X3 < maxC(x||), these forms of 
the total field will not generally be valid. In particular, 
when we are above the surface but still below its max- 
imum point, i.e., C(''^||) ^ < maxC(x||), the expres- 
sion for the scattered field will also have terms contain- 
ing exp [iq|| -xii — iai(gj|)a::3] . Similarly, the transmitted 
field in the surface region has to contain an additional 
term similar to Eq. (10a) but with the exponential func- 
tion replaced by exp [iqy -xy -I- ia2('?|| )2^3] (and a differ- 
ent amplitude). 



If the surface roughness is sufficiently weak, however, 
the asymptotic form of the fields, Eqs. ([s]), Q, and (10 1, 



can be assumed to be a good approximation to the total 
electric field in the surface roughness region. This as- 
sumption is known as the Rayleigh hypothesis [lOlllllITT] . 
in honor of Lord Rayleigh, who first used it in his 
seminal studies of wave scattering from sinusoidal sur- 
faces [ini E]. For a (one-dimensional) sinusoidal sur- 
face, X3 = Cosin(Aa;i), the criterion for the validity of 
the Rayleigh hypothesis, and thus equations that can be 
derived from it (like the reduced Rayleigh equation to be 
introduced below), is known to be (0^ < 0.448, indepen- 
dent of the wavelength of the incident light |3H1 SHI ■ For 
a randomly rough surface, however, the absolute limit of 
validity of this hypothesis is not generally known, though 
some numerical studies have been devoted to finding the 



region of validity for random surfaces |50j . Even if no ab- 
solute criterion for the validity of the Rayleigh hypothe- 
sis for randomly rough surfaces is known, it remains true 
that it is a small-slope hypothesis. In particular, if the 
randomly rough surface is characterized by an rms height 
(5, and a correlation length a (see Sec.[ll|and Ref. [19] for 
details), there seems to be a consensus in the literature on 
the Rayleigh hypothesis being valid if S/a <^ 1 |T71 [5D] . 
We stress that the validity of the Rayleigh hypothesis 
does not require the amplitude of the surface roughness 
to be small, only its slope. 



B. The Reduced Rayleigh Equations 

Under the assumption that the Rayleigh hypothesis 
is valid, the total electric field in the surface region, 
minC(x||) < 2:3 < maxC(x||), can be written in the form 
given by Eqs. ([5]), ([7| and ([lO]) [with Eqs. ^ and 
Hence, these asymptotic fields can be used to satisfy the 
usual boundary conditions on the electromagnetic field 
at the rough surface X3 = C(x||) [SI1I12]- In this way, one 
obtains the so-called Rayleigh equations, a set of coupled 
inhomogeneous integral equations, which the reflection 
and transmission amplitudes should satisfy. 

In the mid-1980s, it was demonstrated by Brown et 
al. [34] that either the reflection or transmission ampli- 
tude could be eliminated from the Rayleigh equations, 
resulting in an integral equation for the remaining am- 
plitude only. Since this latter integral equation con- 
tains only the field above (below) the rough surface, it 
has been termed the reduced Rayleigh equation for re- 
flection (transmission). Subsequently, reduced Rayleigh 
equations for two-dimensional film geometries, i.e., a 
film of finite thickness on top of an infinitely thick sub- 
strate, where only one interface is rough, was derived by 
Soubret et al. [3ll |40] and Leskova [43l|44|. Moreover, 
reduced Rayleigh equations for reflection from clean, 
perfectly conducting, two-dimensional randomly rough 
surfaces [53' and reduced Rayleigh equations for trans- 
mission through clean, penetrable two-dimensional sur- 
faces [32] have been derived. 

For the purposes of the present study, we limit our- 
selves to a scattering system consisting of a clean, 
penetrable, two-dimensional rough surface = C(x||) 
(Sec. [n]). If the scattering amplitudes are organized as 
the 2x2 matrix 



R(qii|k||) = 



i?pp(q|||k,|) i?p,(qj||k||) 
i?^p(q|||k,|) i?s3(qj||k,|) 



(13) 



the reduced Rayleigh equation (for reflection) for this 
geometry can be written in the form 



d^qii /(a2(P||) - ai(q||)|P|| - q||)^,+ , , , / (a2(p||) + ai(k||)|p|| - ky 

,„ " — —. — r —. — \ ^M^(p|| q||)R(q|| kii) = ^ — r 

(27r)2 a2(py)-ai(q||) a2(P||) + ai(ky) 



-M-(p)t|kj,), (14a) 
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where 



and 



= J d^a;|icxp[-i7C(xji)]cxp(-iQji-Xji), (14b) 



M±(-n Pl|9ll±"l('lll)"2(P||)P|rqi| -"«2(P||) [P|| Xq„] \ 

M(P|||q„)-| I (14c) 



where the integrals in Eqs. (14a) and (14b) are over 



the entire qy -plane and xy -plane, respectively. Reduced 
Rayleigh equations for transmission, or film geometries 
with only one rough interface, will have a similar struc- 
ture to Eq. ( 14 ) [SS'.ID] , and can be solved in a completely 
analogous fashion. 

It should be mentioned that the reduced Rayleigh 
equation can serve as a starting point for most, if not 
all, perturbation theoretical approaches to the study 
of scattering from rough surfaces [inj- For example, 
McGurn and Maradudin studied the scattering of light 
from two-dimensional rough surfaces based on the re- 
duced Rayleigh equation, going to fourth order in the 
expansion in the surface profile function, and demon- 
strating the presence of enhanced backscattering j38j . 



C. Mean Differential Reflection Coefficient 

The solution of the reduced Rayleigh equation deter- 
mines the scattering amplitudes i?c(^ (q|| |k|| ). While this 
quantity completely specifies the total field in the region 
above the surface, it is not directly measurable in exper- 
iments. A more useful quantity is the mean differential 
refiection coefficient (DRC), which is defined as the time- 
averaged fraction of the incident power scattered into the 
solid angle dftg about the scattering direction q. The 
mean DRC is defined as 1551 



dns 



1 



cos^ 6*, 



L2 47r2c2 cos 00 



^Q/j(qii|k||)| 



(15) 



herent scattering is 
dRap 



X 



1 



cj^ cos^ ( 



incoh 

|-Ra/3(q|||k||) 



Att^c^ cos Oq 

2 



(i?„^(q„|k||))| 



(16) 



The contribution to the mean DRC from the coher- 
ently scattered light is given by the difference between 
Eqs. (1151) and (fiel). 



D. Conservation of Energy 

As a way to check the accuracy of our results, it is use- 
ful to investigate energy conservation. If we consider a 
metallic substrate with no absorption, the reflected power 
should be equal to the incident power. The fraction of 
the incident light of polarization (3 which is scattered into 
polarization a is given by the integral of the correspond- 
ing mean DRC over the upper hemisphere: 



dRafi 



(17) 



For a non-absorbing metal, if we send in light of polar- 
ization /?, we should have 



^^Q^ — 1, 



(18) 



if energy is conserved. While the conservation of energy 
is useful as a relatively simple test, it is important to 
note that it is a necessary, but not sufficient, condition 
for correct results. 



where is the area of the plane — covered by 
the rough surface. In this work, we are mainly inter- 
ested in diffuse (incoherent) scattering. Since we consider 
weakly rough surfaces, the specular (coherent) scattering 
will dominate, and it will be convenient to separate the 
mean DRC into its coherent and incoherent parts. By 
coherent scattering, we mean the part of the scattered 
light which does not cancel when the ensemble average 
of Rap is taken, i.e., the part where the scattered field 
is in phase between surface realizations. Conversely, the 
incoherent part is the part which cancels in the ensemble 
average. The component of the mean DRC from inco- 



IV. NUMERICAL SOLUTION OF THE 
REDUCED RAYLEIGH EQUATION 

The starting point for the numerical solution of the re- 
duced Rayleigh equation is a discretely sampled surface, 
from which we wish to calculate the refiection. We will 
limit our discussion to quadratic surfaces of size L x L, 
sampled on a quadratic grid of x points with a 
grid constant 



Ax = 



(19) 
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In this paper, we will present results for numerically 
generated random surfaces. These are generated by what 
is known as the Fourier filtering method. Briefly, it con- 
sists of generating uncorrelated random numbers with 
a Gaussian distribution, transforming them to Fourier 
space, filtering them with the square root of the sur- 
face power spectrum .g(k||), and transforming them back 
to real space. The interested reader is referred to, e.g., 
Refs. [231331. 

The next step towards the numerical solution of the 
reduced Rayleigh equation is the evaluation of the inte- 
grals /(7IQ11) defined in Eq. (14b|. These integrals are 
so-called Fourier integrals and care should be taken when 
evaluating them due to the oscillating integrands [54] . 
Using direct numerical integration routines for their eval- 
uation will typically result in inaccurate results. In- 
stead, a (fast) Fourier transform technique with end 
point corrections may be adapted for their evaluation, 
and the details of the method is outlined in Ref. [Mj . 
However, these calculations are time consuming, since 
/(7IQ11) must be evaluated for all values of the arguments 
7 = ai(P||) - a2(q||) and 7 = ai(p||) - a2(k||) [5S]. 

Instead, a computationally more efficient way of evalu- 
ating /(7IQ11) is to assume that the exponential function 
exp [— i7C(x||)] , present in the definition of /(7|Q||), can 
be expanded in powers of the surface profile function, and 
then evaluating the resulting expression term-by-term by 
Fourier transform. This gives 



In practice, the sum in Eq. (20a I will be truncated at a 
finite value n = J, and the Fourier transforms are calcu- 
lated using a fast Fourier transform (FFT) algorithm. 

The advantage of using Eqs. (20) for calculating 
/(7IQ11), rather than the method of Ref. [54], is that 
the Fourier transform of each power of C(x||) can be per- 
formed once, and changing the argument 7 in /(7|Q|[) 
will not require additional Fourier transforms to be eval- 
uated. This results in a significant reduction in com- 
putational time. The same method has previously been 
applied successfully to the numerical solution of the one- 
dimensional reduced Rayleigh equation [35H37] . 

It should be noted that the Taylor expansion used to 
arrive at Eq. (20 1 requires that |7C(xj|)| <C 1 to converge 
reasonably fast, putting additional constraints on the am- 
plitude of the surface roughness which may be more re- 
strictive than those introduced by the Rayleigh hypoth- 
esis. Hence, surfaces exist for which the Rayleigh hy- 
pothesis is satisfied, but the above expansion method will 
not converge, and the more time-consuming approach of 
Ref. [53] will have to be applied. 

Next, we need to truncate and discretize the integral 
over q|| in Eq. ( 14a). We discretize on a grid of equidis- 
tant points, with spacing Aq, such that 



Q 
2 



iAq, 



Q 
2 



(21) 



°° (-i7)" . where i,j = 0, 1,2, ... ,7Vq - 1, and Q = A(7(iV, - 1). 

/(7IQ11) = ^C*'"''(Qll)i (20a) Here, N„ denotes the number of points along each axis 



where C"'(Q||) denotes the Fourier transform of the nth 
power of the profile function, i.e.. 



d2a;||C"(x||)exp(-iQ|| -xn 



(20b) 



of the grid. Additionally, we limit the integration over 
q|| to the region < Q/2. The choice of a circular 
integration domain reduces the computational cost, and 
will be discussed in more detail in Sec. |VI[ Converting the 
integral into a sum by using a two-dimensional version of 
the standard mid-point quadrature scheme, we get the 
equation: 



AqV ^ -^("2(P||) -ai(q||^J|P|| -q||,„) ^ 
^ y ^ ^ ^ ^M+(p|||q,|..)R(q| 

/ qn,<Q/2 ^^^V i-y^Wij' 



/(a2(p||) + ai(k||)|p|| -k|| 



ij 



|k„) = 



(22) 



a2(p|i) + Q!i(k|| 



-M-(p|||k||). 



Here, the sum is to be taken over all qn . . such that we solve separately for either p-polarized incident light. 



q\\i, < Q/2, where gy^ 



This sum yields a ma- 



trix equation where the unknowns are the four compo- 
nents of R(q||^^.|k||). It is evident from Eq. (j8| that if we 
consider incident light of either p or s polarization, we 
need only calculate two of the components of the scatter- 
ing amplitude to fully specify the reflected field. Hence, 



i.e., Rpp and Rsp, or s-polarized incident light, i.e., Rss 
and Rps- In either case, we have twice as many unknowns 
as the number of values of qii . . included in the sum in 



Eq. (22). Note that the left hand side of the equation 



system is the same for both incident polarizations, and 
will also remain the same for all angles of incidence, as 
kii only enters at the right hand side of Eq. (22 1. 
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In order to solve for all unknowns, we need to discretize 
P|| as well, to obtain a closed set of linear equations. Us- 
ing the same grid for py as for qn will give us the neces- 
sary number of equations, as Eq.( 22 ) yields two equations 
for each value of p||. Since we integrate over a circular 
qy domain, with qy discretized on a quadratic grid, the 
exact number of values of q|[^_^. will depend on the par- 
ticular values of Q and iVg, but will be approximately 

In order to take advantage of the method for calculat- 
ing /(7|Qy) described by Eq. (20l, it is essential that all 
possible values of py — qy and p|| — ky [see Eq. (22|] fall 
on the grid of wave vectors Qy resolved by the Fourier 
transform of the surface profile we used in that calcula- 
tion. First, we note that when py and qy are discretized 
on the same quadratic grid, the number of possible val- 
ues for each component of py — qy will always be an odd 
number, 2Nq — 1, where Nq is the number of possible 
values for each component of py and qy. Thus, by choos- 
ing Nq such that 2Nq — 1 equals the number of elements 
along each axis of the EFT of the surface profile we used 



to calculate the integrals in Eq. ( 20b I , we ensure that the 



required number of points is resolved by the FFT [56) . 
Hence, we choose 



Nq^ 



2 



(23) 



where [x\ is the fioor function of x, which is equal to the 
largest integer less than or equal to x. 

Next, we let Aq equal the resolution of the EFT [M] . 
i.e.. 



Aq = 



277 



(24) 



and we let Q be equal to the highest wavenumber resolved 
by the FFT [51], 



Q^Aq[N,/2\. 
In the end, we get the equation 



(25) 



AqV ^{«2(piifc, -ai(qiu,, PiiH-qii„/ ^ 

^ y ^ —r , \ ^M+(pn,,|qi|..)R(q||..|kn ) = 

, V «2(Pi|fc,)-ai(q|u,) "-^UO^h,^ |i„„; 

hi^.N^/^ (26) 



a2(P|Lj + ai(ky_) - y-UO-llmn'^ 



-M-(P||,,|k|| 



kl 



where qy^^-, as well as p 
on the grid given by Eq. (21) 

0,1, 2,..., Nq 



and kii , are defined 

limn' 

with i, j,k,l,m,n = 



1, and where Nq, Aq and Q are given 



by Eqs. (|23j), ([24| and ([25]), respectively. 

Evaluating Eq. (26 1 for all values of py^^ satisfying 

, such that 



P 



kl 



< Q/2, and assuming one value of ky 



^llmri ^ '^Z''' ^^'^ incident polarization j3, results in 
a closed system of linear equations in i?a^(qy^_^. jky^^^^J 
where a = p, s. Repeating the procedure for both inci- 
dent polarizations allows us to obtain all four components 

of R(qiu,mu„)- 

With the reflection amplitudes i?Q^(qy ._^.|ky^^^J avail- 
able, the contribution to the mean differential reflection 
coefficient from the light that has been scattered incoher- 
ently is obtained from Eq. ( 16 1 after averaging over an 



ensemble of surface realizations. 

In passing we note that to avoid loss of numerical pre- 
cision by operating on numbers with widely different or- 
ders of magnitude, we have rescaled all quantities in our 
problem to dimensionless numbers. When considering an 
incoming wave of wavelength A, angular frequency uj, and 
wave vector k, we have chosen to rescale all lengths in our 
problem by multiplying with uj/c, and all wavenumbers 



by multiplying with c/lu, effectively measuring all lengths 
in units of X/2tt, and the magnitude of wave vectors in 
units of Lj/c. 



V. RESULTS 

To demonstrate the use of the formalism for solving the 
reduced Rayleigh equation, the first set of calculations 
we carried out was for two-dimensional randomly rough 
silver surfaces. The surface roughness was characterized 
by an rms height oi S = 0.025A and an isotropic Gaussian 
power spectrum [Eq. ^] of correlation lengths ai — a2 — 
0.25A. In Figs. [2] and [3] we present simulation results 
for the contribution to the mean differential reflection 
coefficients from light of wavelength (in vacuum) A — 
457.9 nm that was scattered incoherently from a rough 
silver surface of size 25A x 25A, discretized into 319 x 319 
points. The dielectric function of silver at this wavelength 
is £2 = — 7.5 + 0.24i, and the angles of incidence were 
6*0 = 18.24° and (j)a = 45°. 

Figure J2] shows the in-plane scattering for this sys- 
tem. The enhanced backscattering peak, a multiple 
scattering phenomenon, is clearly visible, and is as ex- 
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0.008 



(ai?„d/an,,>i.,coh 
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FIG. 2. (Color online) Incoherent part of the mean differ- 
ential reflection coefficient [Eq. (|f6[)], showing only the in- 
plane scattering as a function of outgoing lateral wave vec- 
tor, averaged over 14,200 randomly rough silver surface re- 
alizations. The wavelength (in vacuum) of the incident light 
was A — 457.9 nm, and the dielectric function of silver at 
this wavelength is £2 = —7.5 -I- 0.24i. The surface power 
spectrum was Gaussian [Eq. ([3|], with correlation lengths 
cti = 12 = 0.25A and rms height 5 = 0.025A. The angle of 
incidence was 80 — 18.24°, the surface covered an area L x L, 
where L = 25A, and the surface was discretized on a grid 
of 319 X 319 points. The position of the specular peak (not 
present in the incoherent part) and the enhanced backscat- 
tering peak are indicated by the vertical dashed lines. 



pected strongest in p — p scattering, since p-polarized 
light has a stronger coupling to surface plasmon polari- 
tons [19]. Figure [3] shows the full angular distribution of 
the mean DRC for the same system. In Figs.[3]Ja)-(c) and 
Figs. [3]jd)-(e) the incident light was p- and s-polarized, 
respectively. FiguresJSTc) andJSif) show scattering into 
s-polarization. Figs. pTb) and pTe) show scattering into 
p-polarization and in Figs. [3|^a) and [Sjd) the polariza- 
tion of the scatted light was not recorded. In particular 
from Fig.[3];b), we observe that the enhancement features 
seen in Fig. [2] at angular position 9s ~ —Oq, are indeed 
enhancements in a well-defined direction corresponding 
to that of retro-reflection, and not some intensity ridge 
structure about this direction (as has been seen for other 
scattering systems 45 ). Moreover, the structures of the 
angular distribution of the intensity of the scattered light 
depicted in Fig. [3] are consistent with what was found by 
recent studies by using other numerical methods @5j |46] . 
The results presented in Figs. [2] and [3] were obtained by 
averaging the DRC over an ensemble consisting of 14,200 
surface realizations. 
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FIG. 3. (Color online) Incoherent part of the mean differen- 
tial reffection coefficient [Eq. ( |16[ )], showing the full angular 
distribution as a function of outgoing lateral wave vector. All 
parameters are the same as in Fig. [2] The specular position 
is indicated by the white dots. 



A test of energy conservation was performed by simu- 
lating the scattering of light from a non-absorbing silver 
surface (Im £2 = 0) with otherwise the same parameters 
as those used to obtain the results of Figs. [2] and |3] For 
this scattering system we found — 1| < 0.0003, i.e., 
energy is conserved to within 0.03%, something that tes- 
tifies to the accuracy of the approach and a satisfactory 
discretization. 

As a further test, we studied the scattering from a set 
of (absorbing) silver surfaces with the same parameters 
used to obtain Figs. [2] and [3j except that the rms rough- 
ness S was varied between and 0.045A, while the corre- 
lation lengths were held constant at ai = 02 = 0.25A = a. 
For the purpose of comparison, we also performed sim- 
ulations for a similar set of surfaces but assuming no 
absorption, i.e., we used £2 = —7.5. The results of these 
tests are presented in Fig. |4] 

The reduced Rayleigh equation is only valid for sur- 
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FIG. 4. (Color online) Ratio of reflected power to incident 
power, W, as a function of ratio between rms roughness and 
correlation length, S/a. Surface size and resolution were the 
same as for Fig. [2j and the surface was randomly rough with 
a Gaussian power spectrum, correlation length was kept con- 
stant at a = ai = a2 = 0.25A, while the rms roughness 5 was 
varied from 0.0 to 0.045A. The Fresnel coefficients have been 
included for comparison. 



faces of small slopes [17] . We have found that at least for 
the parameters used in obtaining Fig. |4j our code gives 
good results for an rms roughness to correlation-length 
ratio S/a < 0.12, as judged by energy conservation. For 
larger values of S/a, the results look qualitatively much 
the same, but the ratio of reflected to incident power 
starts to become nonphysical (increasing past 1), as seen 
in Fig. |4] It is noted that decreasing the sampling interval 
Aq, with Q unchanged, did not change this conclusion in 
any significant way, indicating that the observed lack of 
energy conservation was not caused by poor resolution in 
discretizing the integral over q|| . 

The next set of calculations we performed was for a di- 
electric substrate characterized by £2 = 2.64. Otherwise, 
all roughness parameters were the same as for the silver 
surface used to produce Figs. [2] and [3j The mean differen- 
tial reflection coefficient for light scattered incoherently 
by the rough dielectric surface is presented in Fig. |5l By 
comparing these results to those presented in Fig.lSj we 
notice that the dielectric reflects less than the silver (the 
figures show only the incoherent scattering, but the same 
holds for the coherent part), which is as expected. The 
ratio of reflected to incident power for these data was 
U = 0.0467 for p-polarized light at an angle of incidence 
of 0Q = 18.24°. Moreover, from Fig. [5] we also notice the 
absence of the enhanced backscattering peak, which is 
also to be expected since this phenomenon (for a weakly 
rough surface) requires the excitation of surface guided 
modes [19]. Note that for a transparent substrate, it is 
not possible to verify the conservation of energy without 
also calculating the transmitted field. Therefore, energy 
conservation has not been tested for the dielectric sub- 
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FIG. 5. (Color online) The same as in Fig. [s] except that 
£2 = 2.64, and the results are averaged over 21,800 randomly 
rough surfaces. 



strate geometry. 

So far, we have exclusively considered surfaces with 
statistically isotropic roughness. For the results pre- 
sented in Fig. [6| we simulated the light scattering from 
a silver surface of the same parameters as those assumed 
in producing the results of Figs. [2] and |3| except that 
now the surface power spectrum was anisotropic, with 
correlation lengths ai = 0.25A in the xi direction and 
02 = 0.75A in the X2 direction and an rms roughness of 
S = 0.025A. Figure [6] shows the incoherent part of the 
mean DRC averaged over 6,800 surface realizations. In 
this case, there is more diffuse scattering along the xi 
direction than the X2 direction, which is to be expected, 
since a shorter correlation length means the height of 
the surface changes more rapidly when moving along the 
surface in this direction. The interested reader is encour- 
aged to consult Ref . 33J for a more detailed study of light 
scattering from anisotropic surfaces. 

Finally, for the results presented in Fig. [7| we have 
simulated the scattering of light from a surface of size 
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FIG. 6. (Color online) The same as in Fig. |3j except the 
correlation length of the Gaussian roughness, which is ai = 
0.25A in the xi direction and a2 = 0.75A in the X2 direction, 
and the results are the average of an ensemble of 6,800 surface 
realizations. 



25A X 25A, discretized into 319 x 319 points, with £2 = 
— 16 + 1.088i, corresponding to silver at a wavelength 
A — 632.8 nm. The surface power spectrum was cylin- 
drical [see Eq. Q], with fc_ = 0.82a;/c, fc+ = 1.97a;/c 
and rms roughness S = 0.025A, and the angles of inci- 
dence were {9o,(j)o) = (1.6°, 45°). Figure [7] shows the 
in-plane, incoherent part of the mean differential reflec- 
tion coefficient averaged over 7,000 surface realizations. 

From perturbation theory [T7l[T9], we know that for an 
incident wave of lateral wave vector ky to be scattered 
via single scattering into a reflected wave of lateral wave 
vector qy, we must have g{ci\\ — kji) > 0, where ^(ky) is 
the surface power spectrum [Eq. ([2])]. Since the power 
spectrum in this case is zero for |qy — ky| < 0.82w/c, we 
have no contribution from single scattering in the angular 
interval from 9s = —53.5° to 9s = 56.7° (for the angles of 
incidence assumed here). The enhanced backscattering 
peak, which is due to multiple scattering processes, is 
clearly visible in Fig. [t] (at 9s = —Oq) partly because it is 
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FIG. 7. (Color online) Incoherent part of the mean differ- 
ential reflection coefficient [Eq. ( |16[ )], showing only the in- 
plane scattering as function of outgoing lateral wave vec- 
tor, averaged over 7,000 surface realizations with dielectric 
constant £2 ~ —16 + 1.088i, which corresponds to silver at 
A = 632.8 nm. The surface power spectrum was of the cylin- 
drical type [Eq. Q], with fc_ = 0.82a;/c, k+ = 1.97aj/c, and 
rms roughness 5 = 0.025A. The angles of incidence were 
do = 1.6° and Sn = 45°. 



not masked by a strong single scattering contribution. 



VI. DISCUSSION 

A challenge faced when performing a direct numerical 
solution of the reduced Rayleigh equation for the scat- 
tering of light from two-dimensional rough surfaces is the 
numerical complexity. In this section, we discuss some of 
these issues in detail. 



A. Memory Requirements 

Part of the challenge of a purely numerical solution of 
the reduced Rayleigh equation by the formalism intro- 
duced by this study, is that it requires a relatively large 
amount of memory. With approximately JV = (7r/4)iV^ 
possible values for qy , the coefficient matrix of the linear 
equation system will contain approximately (2A/')^ ele- 
ments, where the factor 2 comes from the two outgoing 
polarizations. Hence, the memory required to hold the 
left hand side of the equation system will be approxi- 
mately 4jV^77, where r] is the number of bytes used to 
store one complex number. 
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If each element is a single precision complex number, 
which is what was used to obtain the results presented in 
this paper, then 77 = 8 bytes, and the matrix will require 
approximately bytes of memory for storage. For 

instance, with the choice Nx — 319, which was used in 
all the simulations presented in this paper, and that cor- 
responds to Nq = 160 [Eq. (23l], the coefficient matrix 



will take up approximately 12 GB of memory. 

Note that if we instead performed the qy integration 
in Eq. (14a) over a square domain of edge Q, the number 



of elements in the resulting coefhcient matrix would be 
(2iV2)2 ^ (16/7r2)(27V)2. Hence, by restricting the qy 
integration present in the reduced Rayleigh equation to 
a circular domain of radius Q/2, the memory footprint of 
the simulation is approximately 7r^/16 « 0.62 of what it 
would have been if a square integration domain of edge Q 
was used. For this reason, a circular integration domain 
has been used in obtaining the results presented in this 
paper. However, we have checked and found that using 
a square qy integration domain of a similar size will not 
affect the results in any noticeable way. Indeed, if this 
was not the case, it would be a sign that Q was too small. 

When determining the system size, we can freely 
choose the length of the edge of the square surface, L, 
and the number of sampling points along each direction, 
Nx- These parameters will then fix the resolution of the 
surface, Ax, the resolution in wave vector space, Ag, the 
number of resolved wave vectors, TV^, and the cutoff in 



the q|| integral, Q [see Eqs. ([T9|, ([24]), ([23]), and ([25])]. 
The combination of Aq and Q then determines the num- 
ber of resolved wave vectors that actually fall inside the 
propagating region, |qy| < w/c, which is identical to the 
number of data points used to produce, e.g.. Fig. [3] 

As we are not free to choose all the parameters of the 
system, it is clear that some kind of compromise is nec- 
essary. The number of sampling points of the surface 
along each direction, Nx, and how it determines Nq via 
Eq. (23), determines the amount of memory needed to 



hold the coefficient matrix, as well as the time required 
to solve the corresponding linear set of equations. Hence, 
the parameter Nx is likely limited by practical considera- 
tions, typically by available computer hardware or time. 
For a given value of A^^, , it is possible to choose the edge 
of the square surface, L, to get good surface resolution, 
at the cost of poor resolution in wave vector space, or 
vice versa. Note also that changing L will change Q via 
Ag [see Eqs. (24 1 and (25)]. If Q is not large enough to 
include evanescent surface modes, like surface plasmon 
polaritons, multiple scattering will not be correctly in- 
cluded in the simulations, and the results can therefore 
not be trusted. The optimal compromise between values 
of Nx and L depends on the system under study. 



B. Calculation Time 

The simulations presented in this paper were per- 
formed on shared-memory machines with 24 GB of mem- 



ory and two six-core 2.4 GHz AMD Opteron processors, 
running version 2.6.18 of the Linux operating system. 
The code was parallelized using the MPI library, and 
the setup of the linear set of equations ran on all 12 
cores in the timing examples given. The linear equation 
solver used was a parallel, dense solver based on LU- 
decomposition [54 (PCGESV from ScaLAPACK), which 
runs efficiently on all 12 cores. Setting up the equation 
system scaled almost perfectly to a large number of cores, 
while the solver scaled less well, due to the need for com- 
munication. Numerically solving the reduced Rayleigh 
equation for the scattering amplitudes associated with 
one realization of a rough surface, discretized onto a grid 
of 319 X 319 points, took approximately 17 minutes on 
the architecture described above, and required about 12 
GB of memory. Out of this time, approximately 100 
seconds was spent setting up the equation system, 950 
seconds was spent solving it by LU decomposition, and 
typically around 1 second was spent on other tasks, in- 
cluding writing data to disk. Table [l] shows timing and 
memory details of the calculations, including other sys- 
tem sizes. 

Based on the discussion in Sec. IVI A| we note that the 



use of a circular qy integration domain also significantly 
reduces the time required to solve the resulting linear 
system of equations. When using a dense solver, the time 
to solve the systems scales as the cube of the number 
of unknowns. Thus we expect the CPU time to solve 
the matrix system for a circular integration domain of 
radius Q/2 to be about half {n'^ /2^) the time to solve 
the corresponding system using a square domain of edge 
Q. 

The ratio of the time spent solving one equation sys- 
tem to the total simulation time per surface realization 
increases with increasing system size, as the time to set 
up the equation system is 0{N^), while the time to solve 
the linear system by LU decomposition scales as 0[Nx)- 
It is clear from Table [l] that for any surface of useful size 
the runtime is completely dominated by the time spent 
in solving the linear set of equations. 

Since the time solving the equation system dominates 
the overall simulation time, we investigated if one could 
improve the performance of the simulations by using an 
iterative solver instead of a direct solver based on LU de- 
composition. For example, Simonsen et al. |57j recently 
reported good performance using BiCGStab [S5| on a 
dense matrix system of a similar size. In our prelimi- 
nary investigations into using iterative solvers, we found 
that convergence with BiCGStab was slow and unreliable 
for our linear equation systems. However, it should be 
stressed that we did not use a preconditioning scheme, 
which could potentially yield significantly improved con- 
vergence. 



From Eq. ( 14a ) it follows that changing the angles of 



incidence and/ or the polarization of the incident light 
changes only the right hand side of the equation system 
to be solved. Hence, an advantage of using LU decompo- 
sition (over iterative solvers) is that the additional time 
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TABLE I. Walltime spent to solve the RRE for various values 
of Nx on a shared-memory machine with two six-core 2.4 GHz 
AMD Opteron processors. Included are total time (ttot), time 
to setup the coefficient matrix of the equation system (tLHs) 
and the time to solve the equation system {tsoive)- Also in- 
cluded is the memory required to store the coefficient matrix 
of the linear equation system for each run (Mlhs)- The time 
to set up the right hand side of the linear equation system is 
negligible compared to the left hand side, and have therefore 
not been included here. 





iLHs(s) 


tLu(s) 


ttot(s) 


AfLHs(GB) 


199 


10 


58 


69 


1.8 


239 


28 


171 


200 


3.8 


279 


56 


429 


486 


7.0 


319 


97 


946 


1,045 


12.0 


369 


154 


1,916 


2,074 


19.2 


399 


266 


3,625 


3,895 


29.4 



required to solve the system for several right hand sides 
is negligible, since the overall majority of time is spent 
factorizing the matrix. Conversely, the time spent us- 
ing an iterative solver (like BiCGStab) will scale linearly 
with the number of right hand sides. For these reasons, 
we have chosen to use an LU-based solver. 



C. GPU implementation 

Currently, performing simulations like those presented 
in this paper on a single desktop computer is pro- 
hibitively time consuming due to inadequate floating 
point performance. However, the increasing availability 
of powerful graphics processing units (CPUs) has the po- 
tential to provide computing power comparable to that 
of a powerful parallel machine, but at a fraction of the 
cost. As the most time-consuming step in our simu- 
lations is the LU decomposition of the system matrix 
(see Table |l|, this is where efforts should be made to 
optimize the code. With this in mind, the simulation 
code was adapted to (optionally) employ version 1.0 of 
the MAGMA library [SH] for CPU-based LU decomposi- 
tion. Performance was compared between a regular su- 
percomputing service and a GPGPU (General Purpose 
GPU) tcstbcd. On the regular service, the code was run- 
ning on a single compute node containing two AMD 2.3 
GHz 16-core processors and 32 GB of main memory. On 
the GPGPU testbed, the hardware consisted of a single 
Nvidia Fermi C2050 processor with 3 GB of dedicated 
memory and 32 GB of main system memory. For these 
two computer systems, the initial performance tests indi- 
cated that the LU decomposition took comparable time 
on the two architectures for a system of size Nq = 100 
(the difference was less than 10%). The time using the 
GPGPU testbed included the transfer of the system ma- 
trix to and from the Fermi card and the decomposition of 
the matrix. Even though these results are preliminary, 
it demonstrates that there is a possibility of perform- 



ing simulations like those reported in this study without 
having to resort to costly supercomputing resources. In- 
stead, even with limited financial means, they may be 
performed on single desktop computers with a state-of- 
the-art GPU. 



VII. CONCLUSION 

We have introduced a formalism for performing non- 
perturbative, purely numerical, solutions of the reduced 
Rayleigh equation for the reflection of light from two- 
dimensional penetrable rough surfaces, characterized by 
a complex dielectric function e{uj). 

As an example, we have used this formalism to carry 
out simulations of the scattering of p- or s-polarized light 
from two-dimensional randomly rough dielectric and 
metallic surfaces characterized by isotropic or anisotropic 
Gaussian and cylindrical power spectra. From the 
scattering amplitudes, obtained by solving the reduced 
Rayleigh equation, we calculate the mean differential re- 
flection coefficients, and we calculate the full angular dis- 
tribution of the scattered light, with polarization infor- 
mation. For the scattering of light from weakly rough 
metal surfaces, the mean differential refiection coefficient 
shows a well-defined peak in the retro-reflection direc- 
tion (the enhanced backscattering phenomenon). From 
previous experimental and theoretical work, this is to 
be expected for such scattering systems. Moreover, the 
obtained angular distributions of the intensity of the 
scattered light show the symmetry properties found for 
strongly rough surfaces in recent studies using other sim- 
ulation methods. 

For the purpose of evaluating the accuracy of our sim- 
ulation results, we used the conservation of energy for a 
corresponding non-absorbing scattering system. This is 
a required, but not sufficient, condition for the correct- 
ness of the numerical simulations. By this method, we 
found that within the validity of the reduced Rayleigh 
equation our code produces reliable results, at least for 
the parameters assumed in this study. In particular, for 
a rough non-absorbing metal surface of the parameters 
used in this study, energy was conserved to within 0.03%, 
or better. This testifies to the accuracy of the approach 
and a satisfactory discretization. Moreover, we also per- 
formed simulations of the scattered intensity for systems 
where the rms roughness of the surface was systemati- 
cally increased from zero with the other parameters kept 
unchanged. It was found that energy conservation was 
well satisfied (for the parameters assumed) when the ra- 
tio of rms roughness {5) to correlation length (a), satisfied 
5/a < 0.12. 

We believe that the results of this paper provide an im- 
portant addition to the collection of available methods for 
the numerical simulation of the scattering of light from 
rough surfaces. The developed approach can be applied 
to a wide range of scattering systems, including clean and 
multilayered scattering systems, that are relevant for nu- 
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merous applications. 
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